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Abstract. The density matrix renormalization group (DMRG) method and its 
applications to finite temperatures and two-dimensional systems are reviewed. The 
basic idea of the original DMRG method, which allows precise study of the ground 
state properties and low-energy excitations, is presented for models which include 
long-range interactions. The DMRG scheme is then applied to the diagonalization of 
the quantum transfer matrix for one-dimensional systems, and a reliable algorithm at 
finite temperatures is formulated. Dynamic correlation functions at finite temperatures 
are calculated from the eigenvectors of the quantum transfer matrix with analytical 
continuation to the real frequency axis. An application of the DMRG method to two- 
dimensional quantum systems in a magnetic field is demonstrated and reliable results 
for quantum Hall systems are presented. 



Application of the density matrix renormalization group method 



2 



1. Introduction 

Numerical calculations are now commonly used in various studies in physics. The recent 
remarkable advance in computer technology and the development of efficient algorithms 
provide us with a growing opportunity to resolve behaviors of various interacting systems 
leading us to a new understanding of natural phenomena. In condensed matter physics, 
one of the recent progresses has been made by the density matrix renormalization group 
(DMRG) method, which was developed by Steven White in 1992 PQEj- This method is 
a type of variational method combined with a real space renormalization group method, 
which enables us to obtain the ground-state wavefunction of large-size systems with 
controlled high accuracy. The DMRG method has excellent features compared with 
other standard numerical methods. In contrast to the quantum Monte Carlo method, 
the DMRG method is free from statistical errors and the negative sign problem which 
hampers convergence of physical quantities at low temperatures. Furthermore, the 
essentially exact ground state of large systems extending the limitation of the exact 
diagonalization method is obtained within a restricted number of basis states. The 
truncation error caused by the restriction of basis states is systematically controlled 
by the density matrix calculated from the ground-state wavefunction and the obtained 
results are easily improved by increasing the number of basis states retained in the 
system. 

The ground-state properties of frustrated quantum spin systems and the strongly 
correlated electron systems in one dimension have been extensively studied by the 
DMRG method. Our understanding of these systems was greatly advanced and 
such successful calculations on one-dimensional quantum systems have promoted many 
applications to other problems. In 1995, Nishino applied the DMRG method to the 
transfer matrix and he calculated the partition function of two-dimensional classical 
systems at finite temperatures In 1996, Bursill et al applied this method to 

the quantum transfer matrix (QTM) for S = 1/2 XY spin chains, and calculated 
thermodynamic quantities at finite temperatures. However, they used symmetric 
projection |V aH )(V aR | for the calculation of the density matrix and reliable results were 
obtained only down to T/J = 0.2 Correct calculations were done for S = 1/2 
Heisenberg spin chains by Wang and Xiang 5\, and Shibata [Hj in 1997, and they 
obtained accurate results down to T/J = 0.01. The thermodynamic quantities and spin 
correlation length of the Heisenberg chains with S = 1/2, 1, and 3/2 were systematically 
calculated by Xiang [7j. The stability of this method for fermion systems was shown 
in the calculation of thermodynamic quantities of the one-dimensional Kondo lattice 
model, which is a canonical model of heavy fermion systems [HJ. In 1999 the reliability 
of the method was also shown for frustrated quantum spin chains, t-J ladders and 
the Kondo lattice model away from half-filling [HI E3 II lj . where the negative sign 
problem arises in the quantum Monte Carlo simulations. Dynamic correlation functions 
and one-particle Green's functions in the Kondo lattice model were calculated at finite 
temperatures, and the temperature-induced gap formation in the Kondo insulators was 
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clarified ^2]- There has been a review on the DMRG calculations of thermodynamic 
and dynamic quantities of the Kondo lattice model |13j . Autocorrelations in quantum 
spin chains were systematically studied by Naef et al and its reliability was discussed 
[T3] . The thermodynamic quantities of a spin-Peierls system of CuGeC>3 was studied 
by Kliimper et al . Rommer and Eggert applied the finite-temperature DMRG 
method to impurity Kondo problems, and their results were shown to be consistent 
with the predictions of field theory calculations [TH] . The quasi-particle density of states 
and the dynamic spin and charge correlation functions of the doped one- dimensional 
Kondo insulators were calculated and their doping dependence was clarified |17j . In 
2000, Naef and Wang ^H] studied the nuclear spin relaxation rate 1/T\ in the two- 
leg spin ladder and compared the results to theoretical predictions and experimental 
measures. Maeshima and Okunishi ^H] calculated thermodynamic quantities of 
frustrated quantum spin chains under a magnetic field and studied behaviors near the 
critical fields in the ground-state magnetization process. The magnetization process of 
the Heisenberg spin ladder was analyzed by Wang and Yu (20] , and they determined the 
phase diagram consisting of disordered spin liquid, the Luttinger liquid, spin-polarized 
phases, and a classical regime. Ammon and Imada [2U 1221 EHj systematically studied 
the low-temperature properties of doped £7 = 1 spin chains and clarified behaviors of 
various correlation functions. In 2001, the thermodynamic properties of the £7 = 1/2 
Heisenberg chain with a staggered Dzyaloshinsky-Moriya interaction were studied and 
anomalous behaviors in Yb4As3 were explained [2UI2S]- The thermodynamic properties 
of the t-J chain were studied by Sirker and Kliimper, and the crossover behavior to 
Tomonaga-Luttinger liquid was shown [2E1 121] • Maruyama et al |2H] studied properties 
of a non-magnetic impurity in Kondo insulators. Recently, Sirker and Khaliullin [22] 
studied dimerization in a one-dimensional spin-orbital model with spins £7 = 1. 

The application of the DMRG method to two-dimensional quantum systems is 
currently the most challenging subject and many algorithms have been proposed. 
Most of them use mappings on to effective one- dimensional models with long-range 
interactions, and the standard DMRG method is applied to the effective one-dimensional 
systems. In 1996, White [3U] applied the DMRG method to a two-dimensional frustrated 
quantum spin model for CaV^Og and showed the existence of the spin gap. In 1998, 
White and Scalapino [HI] studied the two-dimensional t-J model with a hole doping and 
found a striped phase. They also studied the competition between stripes and pairing 
in an n-leg t-t'-J model [S2], and made a critical analysis on the phase separation and 
stripe formation in the two-dimensional t-J model [33]. In 2001, Shibata and Yoshioka 
[HUES] applied the DMRG method to quantum Hall systems and the ground state phase 
diagram of two-dimensional electrons in a high Landau level was determined. Xiang et 
al jHH] proposed an efficient mapping on to a one-dimensional system which retains the 
topological characteristics of two-dimensional lattices. They applied this method to the 
£7 = 1/2 Heisenberg model on both square and triangular lattices. In 2003, the ground 
state phase diagram of two-dimensional electrons in the lowest and the second lowest 
Landau levels were determined and the existence of various quantum liquids and charge 
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Figure 1. The system divided into two blocks Bl and Br, whose basis states are 
represented by indices i and j. The solid circles represent local orbitals or quantum 
spins. 

ordered states, such as Laughlin state and the Wigner crystal, was confirmed with new 
stripe states [37|. 

Another approach to higher dimensions was proposed by Xiang in 1996 in which 
the Hamiltonian is renormalized in momentum space [2E1 EH]- This method is clearly 
efficient for weakly interacting systems. The DMRG method has also been applied to 
classical systems in higher dimensions jSj- In 1996 the DMRG method was used to 
renormalize the corner transfer matrix and the corner transfer matrix renormalization 
group (CTMRG) method was formulated [HI]. This method was extended to three- 
dimensional classical models in 1998 [H]. Recently the CTMRG method was extended 
to study one- dimensional stochastic models [12] • The DMRG algorithm has also been 
used to optimize tensor product states in three-dimensional classical systems, and this 
method was applied to two-dimensional quantum spin systems at finite temperatures 
|43j . There has been an analytical study on the spectra of the density matrices used in 
the DMRG calculations for two-dimensional systems |44j . Several topics of the DMRG 
calculations are summarized in a book published in 1998 [4"H] . 

In this topical review, the applications of the DMRG method to finite temperatures 
and two-dimensional quantum systems are reviewed with various techniques in the 
calculation. We first briefly summarize the basic idea of the DMRG method and explain 
how we treat long range interactions. We then apply the DMRG scheme to the QTM and 
formulate an algorithm for the calculations of thermodynamic quantities and dynamic 
correlation functions at finite temperatures. We also apply the DMRG algorithm to two- 
dimensional quantum systems in a magnetic field and calculate various ground states 
realized in quantum Hall systems. 

2. Zero temperature DMRG 

A numerical study of quantum many-body problems requires the handling of 
exponentially increasing basis states of the system. For example, when we represent 
the Hamiltonian of two-spin system with S = 1/2, four basis states, | t, |>, | |, |>, 
| |, |>, and | I, l> are needed. As is easily expected, the dimension of the Hamiltonian 
exponentially increases as 2 N with the increase of the number of spins N. This 
exponential dependence prevents exact numerical diagonalizations for more than 30 
spins in usual cases. The DMRG method is designed to diagonalize the Hamiltonian of 
large-size systems extending this limitation of the exact diagonalization by restricting 
the basis states. Here we briefly review how the DMRG method restricts the basis states 
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2.1. Restriction of basis states 

Let us consider a system consisting of two blocks, Bl and Br (see figureHJ), and represent 
the ground-state wavefunction using the basis states of the blocks described by \i) 
and 

i*>=£*«i«m (i) 

Here we define the following density matrix p for the block Bl, whose basis states are 
represented by |z): 

/4 = E%' (2) 

j 

The norm of the wavefunction is then written by 

(*\*) = Z%*i j = 'E l f& = Tri?. (3) 

i,j i 

The question is how we reduce the number of basis states in the block Bl while keeping 
Tr p L as much as possible. Since Tr p L is equivalent to the sum of the eigenvalues of 
p L , we first solve the eigenvalue equation 

\ , „L „,ct 



W V. 



Y,Pu>v? (4) 



and transform the basis states from \i) to \a) which is defined by 

|a>=£«?IO- (5) 

i 

Since p L is Hermitian, the basis states \a) satisfy the orthogonality 

(«!«') = £«)*<' = <W- (6) 

i 

After this transformation, the density matrix becomes diagonal and the norm of 
the wavefunction is written by 

= Tr p L = (7) 

a 

This result clearly shows that the norm is efficiently preserved by keeping only 
eigenvectors whose eigenvalues w a are large. Thus the optimal m basis states are 
obtained from the eigenstates of the m largest eigenvalues of the density matrix. Since 
p L is the product of and its Hermitian conjugate the eigenvalues w a are all 
positive and we can sort w a as 

w 1 > w 2 > w 3 > ... > w N > (8) 

where N is the number of basis states in Bl. The truncation error, Er(m), of the 
wavefunction represented by m basis states in Bl is then given by 

N m 

Er(m) = £ w a = 1 - £ w a . (9) 

a=m+\ a=l 



Application of the density matrix renormalization group method 



6 




120 160 



a 



Figure 2. Eigenvalues of the density matrix obtained in the one-dimensional Kondo 
lattice model. The number of conduction electrons and localized spins are both 32. 
The electrons and the spins are coupled through on-site antiferromagnetic exchange 
coupling J. t is the hopping integral. 
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Figure 3. Schematic diagram for infinite system algorithm of the DMRG. 



Here we have used (^|^) = J2 a =i wa = !• Thus the truncation error depends on the 
distribution of the eigenvalues of the density matrix, and high accuracy is obtained when 
w a decays rapidly with increasing a. A typical example of w a obtained for the one- 
dimensional Kondo lattice model is shown in figure |2j Even though this model consists 
of conduction electrons and localized spins, the accuracy of 1CT 7 is attained by keeping 
only 150 states in each block. Accurate results in the DMRG calculations are obtained 
when w a decays almost exponentially as shown in figure El 

2.2. Infinite system algorithm 

Here we describe how we obtain the Hamiltonian of a large system within a restricted 
number of basis states. In the infinite system algorithm of the DMRG method[lJ, we 
start from a small system and iteratively extend the system by adding new sites between 
the two blocks, as schematically shown in figure El At each extension of the blocks we 
restrict the basis states using the density matrix defined in section 2.1. In the following, 
we describe the algorithm of the DMRG method applied to one-dimensional quantum 
systems which have long-range interactions. 

We usually start the calculation from a system consisting of four single sites with 
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open boundary conditions. The left and right blocks contain the single site at one end 
of the system, as shown in figure EJ The Hamiltonian is written by 

H = H Bl + H Cl + H° 2 + H Br 

_i_ ^jB L —Ci _|_ j^Cx—Ci _|_ jjC2-B r 

+ h Bl ~ C2 + h Ci ~ Br + h Bl - Br (10) 

where H Bl and H Br contains operators represented by only the basis states of the left 
and right blocks, respectively, and H Cl and H° 2 are the single-site Hamiltonian between 
the two blocks. The interaction and hopping terms between different blocks and single 
sites are represented by H Bl ~ Ci + H Cl '° 2 + H° 2 ' Br and H Bl ~ C2 + H Ci ~ Br + H Bl ~ Br , 
where the former three terms describe the nearest-neighbor interactions and hoppings, 
and the latter three terms represent next-nearest-neighbor and further long-range 
interactions and hoppings. In order to avoid complexity we have assumed that each 
inter-site interaction and hopping term is written by operators at two different sites, and 
there are no multi-site interactions such as Sr (S m x S„) with I ^ m ^ n. The Heisenberg 
model with nearest-neighbor interaction J\ and next-nearest-neighbor interaction J2 in 
magnetic field h z is written by H Bl = h z Sf, H Cl = h z S z 2 , H c ' 2 = h z SI, H Br = h z S(, 

H B L - Cl = JiSi . H C X -C 2 = JiS2 . R C 2 -B R = JiSs . R B L -C 2 = JsSi . 

H d-B R = J 2 g 2 . g 4j and H B L -B R = q 

We then calculate the ground state- wavefunction \^) by diagonalizing the above 
Hamiltonian. The obtained wavefunction is represented in the basis states of the blocks 
and single sites as 

l*>= E ^^ L ^^^M\^l)\^2)\^R) (H) 

where and are the basis states of B L and B R , and and \i 2 ) are the basis 
states of C\ and Ci. We next extend the blocks by including the neighboring single 
site. The basis states of the extended blocks Bl ® C\ and Br ® C2 are represented 
by the product space of the original block and the single site, |*z,)|«i) and Izr)^)- We 
then restrict the basis states of the extended blocks. For this purpose we calculate the 
density matrices for the extended blocks, p L and p R defined in section 2.1 

Pi L ii,i' L ^ = E ^iLiii2iR&i' L i[i2i R ( 12 ) 



Pi R i2,i' R i' 2 E ^ ihilWR^ i L iii' 2 i' R - (13) 

We diagonalize the density matrices to obtain the eigenvectors v a and the corresponding 
eigenvalues w a which satisfy 

™ aL <\ = E4, iW ^; ( 14 ) 
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The new basis states of the extended blocks are obtained from the eigenvectors of the 
m largest eigenvalues \oll) and |or) defined by 

M= E<nl^)Ki) (16) 

l«*> = ECJ^)- (17) 

We then represent all operators using these new basis states. The spin operator at the 
first site (S()i Li ' L in the left block Bl is given by 

(s?)«w= e <\(<imki £ '- (is) 

Similarly, the spin operators at the second site which has been included in the 

new block (Bl) new IS represented by 

(s z 2 )* LaL >= E <\(<iO*(sf)w- (is) 

The block Hamiltonian (H BL ) new = H Bl + if * 1 + H Bl ~ Ci and the interaction between 
the new block and the single site (H L 1 ) new = H 1 2 + H Bl ~ c ' 2 are obtained by 

(20) 

Similarly, (H Ba ) new and (H C2 ~ BR ) new are given by 

( H aR a > R )new = E ^fllaKvia')* { H fj> R 8i 2 i' 2 + ^\^i R i' R + H ° R %fj>J 

(22) 

(^Sii^i'^aew = E ^laKvia')* {^fX^^R^R + ^Si^t^J ' ^ 

The long range interactions (H BL ~° 2 ) new , (H Cl ~ BR ) ncw and (H BL ~ BR ) new are not simply 
obtained by transforming the previous Hamiltonian, because they are defined after we 
add new single sites between the two blocks. These terms are constructed from the 
operators represented by the basis states of the new blocks and the single sites added 
in the system. The new Hamiltonian of six sites is then written by 

(H) new = (H Bl ) qcw + H c ^ + H° 2 + (H Br ) qbw 

+ (H L - C1 ) ncw + H c ^° 2 + (H C2 ~ BR ) new 

+ (#^- C2 ) new + (# Cl - B *)new + (H B ^ BR ) ncw . (24) 

This Hamiltonian has similar structure to the previous one written in equation (fT0|) . 
and we repeat the above procedure until we obtain the desired size of system. Since 
there is no limitation on the size of system, this algorithm is called the infinite system 
algorithm of the DMRG. 
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Figure 4. Schematic diagram for finite system algorithm of the DMRG. 



2.3. Finite system algorithm 

In the infinite system algorithm of the DMRG we calculate the Hamiltonian of any 
size of systems by iteratively extending the blocks. However, in order to deal with 
exponential increase in the number the basis states, we truncate the basis states. This 
truncation introduces an artificial error in the wavefunction, and we need to minimize 
the error by reconstructing the basis states within a fixed size of system. The algorithm 
to find the optimal basis states in a finite system is called the finite system algorithm 
of the DMRG. Here, we describe how we improve the ground-state wavefunction using 
the finite system algorithm of the DMRG. 

Suppose we have extended the system up to 2N + 2 sites by using infinite system 
algorithm of the DMRG. The left and right blocks now contain N sites and the system 
is written by 

B L (N).»B R {N). (25) 

Here Bl(N) and Br(N) represent the left and the right blocks which contain N sites 
and • represents a single site between the two blocks. In order to reconstruct the right 
block first, we extend only the left block and reduce the right block, as shown in figure 
IU We calculate the optimal basis states for Bl(N + 1) by diagonalizing the density 
matrix calculated from the ground-state wavefunction of the system Bl(N) • •Br(N), 
and we obtain the Hamiltonian H Bl for the left block Bl(N+ 1) by a similar calculation 
to that in the infinite system algorithm of the DMRG. We do not need to calculate the 
basis states and the Hamiltonian H Br for Br[N — 1) because they have already been 
obtained in the infinite system algorithm of the DMRG. The system of 2N + 2 sites is 
now represented by 

B L (iV + l)..B R (iV-l). 

We repeat the extension of the left block until the right block becomes single site. 
Then, we extend the right block from the single site to reconstruct the basis states in 
Br. Now the basis states of the right block are always determined from the ground- 
state wavefunction of 2N + 2 sites. This is different from the calculation in the infinite 
system algorithm, where the size of system increases with the extension of the blocks. 
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This difference improves the ground-state wavefunction obtained in the infinite system 
algorithm of the DMRG. 

The extension of the right block is continued until the left block is reduced down 
to the single site. We then extend the left block to reconstruct the basis states in B L . 
We repeat such sweeps until the ground-state energy and its wavefunction converge. 
Since we have reduced Hilbert space of the Hamiltonian by restricting the basis states, 
the obtained ground-state energy is always higher than the exact one. The increase 
in the number of the basis states systematically lowers the ground-state energy and 
its value approaches to the exact one. Thus the finite system algorithm of the DMRG 
method is equivalent to a variational method in which a trial wavefunction is constructed 
automatically under the restriction on the number of basis states. 

2.4- Evaluation of physical quantities 

We evaluate various physical quantities from the wavefunction obtained by the DMRG 
calculations. Since the basis states of the wavefunction alter after the reconstructions 
of the blocks, we need the operators represented in the basis states of the wavefunction. 
These operators are obtained by successive transformations of basis states starting from 
the operators defined in the initial single site. At each reconstruction of the blocks the 
operators such as Sf is transformed as follows. If the ith site is the single site next to 
the block B L , (Sf)^' is transformed by 

(s?) aLaL '= E <\Ki'T( s i)^' (26) 

ihhh' 

where extended new left block include the ith. site and has basis states \ai), and vf^ 
is the eigenvectors of the density matrix defined in equation JHJ). After the ith site is 
included in the block B L , S? is transformed by 

(S?) aL « L i= E <\«X)*(^)^'- (27) 

The product of the two operators at sites i and j, such as S*Sj are transformed by 

{S z t S*) aLaL ,= £ vl\{v^ h ,)\St) ihi AS^ (28) 

where the jth site is the single site next to the original left block and the zth site 
is in the left block Bl. If both the ith and jth sites are in the left block Bl, they are 
transformed by 

(S?S;) aLaL ,= £ v?X{vt\T{SlSI) iLiL , (29) 

After we have obtained the operators represented by the basis states of the 
wavefunction ^i L i 1 i 2 i R \iL)\h)\'i2)\'i'R)i we evaluate physical quantities. The local 
quantities {^\Sf \^) and correlation functions (^\S- Sj\^f) are given by 

OWs?!*) = E Kw^^k^^ (3i) 
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Figure 5. Long range hopping between the two blocks. 



where both Si and Sj are assumed to be in the left block. If Sj is in the right block 
while Si is in the left block, the correlation function (^f\S^ Sj\^) is given by 

(*|5?s;i*> = E K'iMA^'iAs^i^ir.^. m 

iLi^'iiWRiR 1 

2.5. Fermion sign 

When an electron hops from one site to another, the fermion sign appears in the matrix 
element between the two basis states. The fermion sign depends on the number of 
electrons between the two sites, and we need to count the electrons. However, the 
number of electrons at each site in the basis states of the blocks is neither one nor zero 
after we have transformed the basis states using the eigenvectors of the density matrix. 
In order to correctly deal with the fermion sign, it is useful to define new creation 
(annihilation) operators, 4 and (c ia ), which include the fermion sign generated by the 
hopping from the ith site to one end of the block. If we align the up-spin electrons and 
down-spin electrons separately, 4 in the left block are defined by 

4 = 4(_i)JW+JW+-+^ (33) 

where L is the number of sites in the left block (see figure EJ) and are the number 
operators of the electron with spin a at the ith site. The operators cj CT in the right block 
are defined by 

c) a = C } a (-1) N ^ +N ^+- +N ^ . (34) 

The matrix element {A;|4 c jvl0 generated by a long-range hopping between the two 
blocks is then given by 

(k\clc ja \l) = (4)^^(^)^H(-l) (JVi+lCT)il+(iVi+2CT)i2 ^^ fc ^ 2 (35) 
where |/) = |/i)|ii)|i 2 )|^)- The coefficient (-l)( N L+^K+( N L+^)h is the fermion sign 
coming from the electrons with spin a in the two single sites between the blocks. In 
this expression we need not calculate the fermion sign coming from the electrons in 
the blocks, and the matrix elements of long-range hoppings between the two blocks are 
obtained by only counting electrons in the two single sites between the two blocks. 

The operators 4 are obtained by successive transformations starting from the 
original operators defined in the single site. When the left block is enlarged by including 
the iith site, 4 i n the left block is transformed as 

(4)cw= E c^nikvC-i)^ (36) 
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where (—1)^1" is the fermion sign generated by the electron at the iith site, and c ia = c ia 
when the ith site is the right end of the left block, i = L. Since the basis states of 
the single site are the eigenstates of the number operator N^, is either one or zero 
depending on the index l^. 

2.6. Conserved quantities 

In the DMRG calculation, most of the computing time is consumed for the 
diagonalization of the Hamiltonian. Thus, it is important to use an efficient method 
for the diagonalization. Here we consider the basis states classified by the conserved 
quantities of the Hamiltonian. By using these basis states we need only diagonalize the 
Hamiltonian within a subspace specified by a set of quantum numbers. For example, 
Heisenberg model without transverse magnetic field conserves the z-component of the 
total spin, and all the eigenstates of the Hamiltonian are classified by total S z . The 
blocks Bl(2) and Br{2) containing two S — 1/2 spins have basis states | f, T>, | T, !>, 
| I, |>, and | I, l> with quantum numbers of S z = 1, 0, 0, and -1, respectively, and the 
basis states of the total system B L (2)»»B R (2) are constructed from various combinations 
of local basis states which satisfy 



When we calculate the ground state of S% ot = 0, there are 10-configurations which are 



represented by (S* L , S Z C1 , S Z C2 , S Z R ) = (1,1/2,-1/2,-1), (1,-1/2,1/2,-1), (1,-1/2,-1/2,0), (0,- 
1/2,-1/2,1), (0,1/2,-1/2,0), (0,-1/2,1/2,0), (-1,1/2,1/2,0), (0,1/2,1/2,-1), (-1,1/2,-1/2,1), 



(-1,-1/2,1/2,1). The number of total basis states is 20, which is lower than 2 6 = 64 of 
full Hilbert space. Since the Hamiltonian H Bl and H Br conserve total S z , they are 
also block diagonal, and H Bl + H Cl + H° 2 + H Br have matrix elements only within 
each subspace specified by (Sg L , Sq v S^ 2 , S z Br ). The matrix elements between different 
subspaces are generated only by the terms H B - Cl + H Cl ~ C2 + H C2 ~ Br + H Bl ~° 2 + 

The density matrix p iV = Y.j obtained from the ground-state wavefunction 

is also block diagonal, because total S z of the basis states j is uniquely determined by 
the basis states % of the density matrix. This means that the density matrix does not 
have any matrix elements between different subspaces of total S z . Since the density 
matrix is block diagonal, the unitary transformation of the basis states is defined only 
within the same subspace of total S z , and the basis states of the new blocks are also 
classified by total S z . 

By using the basis states classified by the conserved quantities such as the z- 
component of the total spin and the number of electrons, the dimension of the 
Hamiltonian, operators, and density matrix are reduced. The reduction of the 
dimensions improves the efficiency and accuracy of the numerical calculation, and 
reduces the memory space needed in the calculation. 




$L + &C1 + $C2 + Sr — 



const. 



(37) 
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Figure 6. Arrangement of the blocks and the single sites for the periodic boundary 
conditions. 
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Figure 7. Eigenvalues of the density matrix for S=l/2 Heisenberg model with 40 
spins under the periodic boundary conditions (PBC) and the open boundary conditions 
(OBC). 



2.7. Periodic boundary conditions 

In the DMRG calculations, the open boundary conditions are usually used. However, the 
existence of the boundaries sometimes causes large boundary effects on the ground state. 
These boundary effects are removed when we impose the periodic boundary conditions. 
In this case, we usually align the blocks and single sites symmetrically as shown in 
figure El The truncation error under the periodic boundary conditions is generally large 
compared with that obtained in the open boundary conditions. To keep the accuracy 
of the results, we need much more basis states in each block. This is actually shown 
in figure [7J which shows slower decay of the eigenvalues of the density matrix for the 
periodic boundary conditions. 

Under the periodic boundary conditions, the left and right blocks are connected 
with each other through both ends of the blocks. The reduction of the accuracy is 
caused by the fact that the interactions around the boundaries between the blocks and 
the single sites such as H Bl ~ Ci and H c ' 2 ~ Br are represented in the product space of 
the block and single site, and they are not renormalized to reproduce the ground-state 
wavefunction within a small number of basis states. When the interactions across the 
boundaries are weak, the accuracy of the DMRG calculation will be improved. 
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3. Finite temperature DMRG 

The zero-temperature DMRG method described in the above section enables us to 
calculate the ground state wavefunction of large-size of systems with controlled accuracy. 
However, in order to study the properties of the bulk system, we need to extend the 
system in real space direction and analyze the asymptotic behaviors. In the finite- 
temperature DMRG method, the thermodynamic quantities of an infinitely long system 
are directly obtained from the eigenvalues of QTM, and temperature dependence is 
obtained in the process of iterative extensions of the QTM in the /5-direction. In 
this section we describe how we calculate thermodynamic quantities using the DMRG 
method applied to the QTM. 



3.1. Transfer matrix method 

The thermodynamic quantities are calculated from the partition function of the system. 
In one-dimensional quantum systems the partition function is obtained from the 
maximum eigenvalue of the QTM [461 HZj- Here we first briefly summarize the way 
to calculate the QTM and the partition function. 

In the following we suppose that the Hamiltonian consists of only on-site and 
nearest-neighbor interactions so that we can divide the Hamiltonian into two parts 
#odd = En=o _1 h 2n +i,2n+2 and H even = En=o -1 h 2n ,2n+i- Here h iti+1 contains on-site and 
nearest-neighbor interactions between the ith and i+lth sites and each hi^+i commutes 
with all other h i+2n ,i + i+2n- 

[^2n,2n+l> ^2n',2n'+l] = [^2n+l,2n+2; ^2n'+l,2n'+2] = 0. (38) 

The partition function is then decomposed into a product of e~^ hi ' i+1 ' M 
Z = Tr e-? H 



as 



lim Tr 

Af-+oo 



lim Tr 

AT -+oo 



e -f3H oAd /M e -/3H eve 



n/M 



L/2-1 

n e 

n=0 



L/2-1 

Ph2n+l,2n+2/M J~J g- 
n=0 



■Ph 2 



i/M 



M 



(39) 
(40) 

(41) 



where M is the Trotter number j4<2 El E0J- 

Here we introduce imaginary time Tj, whose index j takes only integers up to 2M, 
and we represent the partition function in a tensor product of e _/3/li i+1//M under the 
complete basis set for Tj as shown in figure IBf a), where e~P hi < i+1 ' M is represented by a 
hatched square shown in figure E^b). Each e~P hi < i+1 ' M is then represented by the basis 
states specified by i,i + 1 and Tj,Tj+i. We explicitly write the imaginary time indices 
Tj and Tj+i in e _/3hi 



_1//M , and define the tensor e^ h ^^\j M represented by 



E 



t j+i 



E E E 

+ 1 ' T j + l CT i.r,- a i + l.r 



rj+i 



-l,r. 



)((T-cr- L -,\e- f}hi ' i + l/M \cr'cr' )(rr' a' 
- J+1 / \o %o 1+1 1 e \° i° i+il \° i,Ti° i 



(42) 
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Figure 8. (a) Partition function Z for M = 3 represented by a tensor product of 
e -/3oh eve n anc j e -f} h O dd^ where /i ewe n = h 2n ,2n+i and /i odd = /i2n+i,2n+2- (b) Graphical 
representation of e, f 3 ° h '- z + 1 with do = 8/M. 



where a^ T . is the local variable which specifies the state at the ith site and at imaginary 



time Tj. 



The QTM T n is written by 

A/-1 

T n = lim TT ( e -^n+l, 2 n+2/M -fihwn+i/M. ^ 
M->00 11 (^2j+2,T2j+l) (T2j+l,T2j) > V / 

as a tensor product of e^^^ M with the periodic boundary conditions in the (3- 
direction, t 2 m = t . This matrix T n is graphically represented in figure El When 
the system is translationally symmetric, we can omit the site index n in T n , and the 
partition function is given by 

Z = Tr T L/2 . (44) 

Using the eigenvectors of the transfer matrix as basis states, we can represent the 
partition function using their eigenvalues Aj as 

Z=£Af /2 . (45) 

i 

In the thermodynamic limit L — > oo, Z is determined only by the maximum eigenvalue 

^max • 

Z = \ L ±- (46) 
The free energy par site is then calculated from the maximum eigenvalue X ma x as 

F= - In Z/(L/3) (47) 
= -\n\ max /(2f3). (48) 

The thermodynamic quantities such as the entropy, specific heat and the magnetic 
susceptibility are given by taking the derivative of F with respect to the temperature 
and magnetic field. 
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Figure 9. Quantum transfer matrix T represented by a tensor product of e P h i,*+i/ M 
with M = 3. 



Application of the DMRG to the quantum transfer matrix 

In order to calculate low temperature properties, we need to increase the Trotter 
number M keeping (3/M = At sufficiently small. However, the dimension of the QTM 
exponentially increases with increasing the Trotter number. In this section we describe 
how we apply the DMRG method to QTM and obtain the maximum eigenvalue for large 
M. 

In the zero-temperature DMRG method, we extend the Hamiltonian in the real 
space direction by restricting the number of basis states by using the eigenvectors of the 
density matrix calculated from the ground-state wavefunction. In the finite-temperature 
DMRG method, we extend the transfer matrix in the /^-direction by restricting the 
basis states using the density matrix calculated from the eigenvector of the maximum 
eigenvalue \ max . 

We start from a small QTM with M = 2 and divide it into two blocks, 7" A ( M=1 ) 
and T B( - M=l K as 



j-A(M=l) 

(T2,Tl,To) 



rj-B(M=\) 

(t2,ti,t ) 

— /3o/»odd „-Poh ev , 
(t2,ti) (ti.tq) 



(49) 



where /i dd 



^2n+i,2n+2 and h 

q-(2M) 



h,2n,2n+i- The transfer matrix T is defined by 

(50) 



j-B(M) <j-A(M) 

('t)i T (4M-l,---,2M+2,2A/+l)!' r 2M ) ( T 2M,T(2M-1 



,2,1)> 7 



J-(2M+1) _ q- 



B(M+l/2) 



A(M+X/2) 



■ _ ! (51) 

"(^0,T'(4M+l,---,2A/+3,2M+2)>' r 2Af+l) (t"2M+1 ,T(2M,- •■,2,1) > T o) ^ ' 

depending on the parity of total Trotter number, and we have imposed the periodic 
boundary conditions, ov 4M = o"i )T0 , in the /^-direction. We extend T A and T B as 



-A(M) 

'(T2M+1 ,T2m) (V2M ,T(2M-1, 

-B(M) 



■,2,l), r o) 
Aj^odd 



T ' e 

( T 2M+2;' r (2M+l,---,4,3) i r 2 



(T2Af+l,T(2A/,...,2,l)i' r o) 

j-B(M+l/2) 

(t2M+2,T(2M+1,---,3,2)) t i) 



( , j- A(M +1/2) 

( r 2M+2. T 2A/+l) ( T 2A-/+1 ; T (2Af,---,2,l) i T ()) 



"Aj^odd 
T~2M+ 

-B(M+l/2) 



X 



A(M+1) 



T - < 

(T2M+2,T(2M+2,...,3,2)in) (TL,To) 



(T2M+2;T(2M+1 

j-B(Af+l) 



2,l)! r o) 

(T"2Af+2,T(2Af+l,-.-,2,l))' T "o) 



(52) 
(53) 
(54) 
(55) 
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Figure 10. Schematic diagram for the expansion of the QTM. 



as shown in figure Ell 

At each expansion of the QTM in the /3-direction, we restrict the basis states of 
T A and T B while keeping the accuracy of A max and corresponding eigenvectors, which 
satisfy the following equations 

V 3 L T jf = \ max Vf (56) 



E r^v} 



33 
R 



(57) 



where Tjy represents the matrix element of T between the basis states j and j', and 
V L and V R are the left and the right eigenvectors of the maximum eigenvalue A max , 
normalized by (V L \V R ) = 1. V L and V R are usually different because T is non- 
Hermit ian. 

Since the maximum eigenvalue \ max of T is equivalent to (V L \T\V R ), we need to 
keep the norm (V^V^) as much as possible when we restrict the basis states. To find 
the optimal basis states of T A , we calculate the following density matrix defined by 



Pkk> 



(5? 



where index k represents the basis states correspond to the imaginary time between 
T\ and T2M (between n and t^m-i) for T A ( AI+1 / 2 ^ (T A ^), which is transformed to 
new basis states, and I represents the basis states for the remaining imaginary time: 
{j} = {k, I}. Since (V^V^) is given by the sum of the eigenvalues w a of the density 
matrix, (V L \V R ) = Tr p A = J2kPkk = J2a w a, the optimal basis set within a fixed 
number m of basis states are obtained from the eigenvectors of the m largest eigenvalues. 
The new basis states are thus defined by the eigenvectors v a ( L > and v a ( RS ) of the density 
matrix, which satisfy the eigenvalue equations 



£ 



V 



Pkk' 



w a v 



a(L) 



(59) 
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Figure 11. Eigenvalues of the density matrix p calculated from the eigenvectors for 
Amaz of the QTM in one-dimensional Kondo lattice model at half-filling. 



Pkk' 



.MR) 



a a(R) 
W V k ■ 



(60) 



By using v a ^ and v a ( R *> we transform T A as 



aa' 



E 

kk> 



a {R) T A a'(L) 
I kk' v k' 



(61) 



Here, we have omitted the indices of the basis states corresponding to the part of the 
Hilbert space with imaginary time r and t 2 m+i ( r o an d t 2 m) for 7~ A ( A/+1 / 2 ) (T A ( M ^), 
which are included in T A but independent of the transformation. As a result of the 
non-Hermitian density matrix, the left and the right eigenvectors, v a ^ and v a<yR \ are 
different, and this difference makes a different transformation between the left and right 
basis states of the transfer matrix. Since the left and right eigenvectors v a ^ and v 01 ^ 
are dual orthogonal 



E 



h Uh. 



(62) 



the orthogonality of the left and right basis states is retained. 

In this way, we restrict the number of basis states while keeping the accuracy of 
the maximum eigenvalue and its eigenvectors. The truncation error in the calculation 
is estimated from the sum of the eigenvalues of the density matrix which are truncated 
off, and it is given by 1 — where w k is the kth largest eigenvalue. In figure HP 

we show a typical example of the eigenvalues of the density matrix. We can see only 
a small number of eigenstates have large eigenvalues and the truncation error is about 
10 -4 for m = 50. In such a case, we can repeat the expansion of the QTM and obtain 
its maximum eigenvalue with controlled accuracy. 



3.3. Quantum transfer matrix 



Here we explicitly write the QTM. The QTM is obtained by extracting the matrix 
related to the site indices 2n, 2n + 1, and In + 2 from the partition function decomposed 
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into a tensor product of e, /3 ° /l2 "+ 1 - 2 ™+ 2 and e, /3 ° /l2 "- 2 "+ 1 . We describe these local tensors 

1 (T2m+2,T2m+l) (T2m+l,T2m) 

using the local variables <7j jTj at site i and imaginary time r 3 -. a represents up-spin state 
or down-spin state for S — 1/2 Heisenberg model. The transfer matrix T^ M ^ is then 
written by 

j-(M) = ^ ©(<Ti 1 7 t ,<T2,7 l )ki,7a M _i<7 , 2,7a M _i) • • • 

0-l, T - () ---0-l,T 2M _ 1 

)0( ) (63) 

where 

©(^,^■+1^+1,^+1 1 o- i:Tj a i+ljT .) = (ai,r 3+1 ,^i+i,T 3+1 |e (r *^YI ai ^' <Ji + 1 ^) ( 64 ) 
and fa = P/M = 2(r j+1 - tj) = At. 

The matrix elements Q(o~i jTj+l o~i+i jTj+l \o~i :Tj o~i + i :Tj ) are determined by using 
eigenstates of the local Hamiltonian /i^+i. In the case of 5=1/2 Heisenberg model, 
hij+i is JSj'Sj+i. The eigenvalues E m are (— 3J/4, J/4, J/4, J/4) and the corresponding 
eigenvectors are 



\Vi) = 

\v 2 ) = 

\Vs) = 
\V 4 ) = 



^(IU>-|U» 



V2 

IU> 
I!,!). 



(I U> + l U» 



(65) 

(66) 

(67) 
(68) 



The matrix element (fi, U+i |e ^m+i] | i; is then given by 



_( e /3o3J/4 + e -A)J/4^ 

2 



(69) 



Similarly, off-diagonal elements and other diagonal elements are obtained as follows: 



(U,h + i\e- (3ohi > i+1 \h,h + i) = (UU+i \e 



±( e 0o3J/4 _ -A.J/4N 

2 V ' 

A)^i.i+i| I . I 
I +n 4- 



= e 



-A J/4 



Thus, the matrix 0(<7i jTj+1 (7i + i jTj+1 |<7j jTj (7j + i )T .) is represented by 

/ e" A ' J / 2 

cosh(-/? J/2) sinh(-/5 J/2) 

sinh(-/5 J/2) cosh(-/3 J/2) 

v e -M2 

with the basis states of the matrix 



(70) 
(71) 



(72) 



( I TijTj) Ti+l,Tj)j | TijTj) -U+1,Tj)j I ii,Tj 1 Ti+l,Tj) J I ii,Tjl il+ljTj) )■ 



,T 2/ 
!>T 2 
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In order to explicitly represent the QTM, we need to use the basis states 

( I Ti,Tj-j Ti,Tj-+i/> I Ti,Tj; li,"rj-|_i)) | ■iijTj-) Ti,Tj + i); | ti,Tj j li,Tj--)-i) ) 

defined along the /3-axis at a given site i. The matrix defined in these basis states is 
related to the matrix as 

0(cri,r i C"j,r J+1 |o"j+l,r J O"i+l,r 3+1 ) = ®{&i,T :i+1 (Ti+l,Tj +1 \o'i,TjO'i+l,T j )- (73) 

Thus, the matrix is represented by 

/ e~^ J / 2 cosh(-A,J/2) \ 

sinh(-/3 J/2) 

sinh(-/3 J/2) 

V cosh(-/5 J/2) e-^ J l 2 j 

This is the smallest unit of the QTM shown in figure IHfb). 

The initial blocks of the transfer matrix 7" j4 ( M=1 ) anc l r fB{M=l) j n e q Ua tion (|49j) are 
written as 

T A(M_1) (er 0)T0 , cr 0)Tl | (7i,to , a 1)T2 \ a 2)Tl , cr 2 , 

— T' B( " M ~ 1 \o- 0tTo ,ao, Tl \&1,t , Cl.ral^.ri, &2, 

,Ti°~l,T2 | cr 2,Ti cr 2,r 2 

). (75) 

< 7 1,T 1 

The QTM T^ M=1 ^ is obtained from 7" j4 ( M=1 ) with the periodic boundary conditions in 
the /3-direction a itT2 = c itT0 : 

T (M_1) (o- , To ,a 0(Tl | cr 2 , To , cr 2lTl ) 

2 

= E E E TA(M=1 H^0,r ,^0,r 1 1^ l,r ,^ 1,^1^2,^,0-2,^ ) n^.r 2 ^, T0 - (76) 

CT 1,T < J l,r 2 CT 2,T 2 1=1 

5.^. Renormalization of the quantum transfer matrix 

The QTM at low temperature is obtained by extending the QTM in the /3-direction. 
By using T A ( M=l ) anc l >j~b(m=i) defined j n equation (|75|) . we first calculate the QTM of 
M = 2. This QTM corresponds to a high temperature of T = l/(Mf3 ) = l/(2/3 ) > 1 
since we need to satisfy /3 = At <C 1. We therefore decrease the temperature by 
iteratively increasing Trotter number M as T = l/(M/3 ) — > T" = 1/((M + l)/? ) 
keeping /3q fixed. In the following we describe detailed calculation for the extension of 
the QTM in the /3-direction. 

As is shown in figure EH the T( M=2 ) is written by T A ( M=1 ) and T B ( M=1 ) as 

j-(M-2) ^j-g^p. Q ^ ? ^ Q ^ ? Q ^ | a 2 ,T ^ (T2 t Ti , <?2,t 2 i a 2,T 3 ) 

= E '^ B( ' M ~ 1 \ (T 0,T2, ffO,T 3 |o"l,T 2 , O"l,T |ff2,T 3 , O"2,T ) 

CT 1,TQ !°"1,T 2 

x 7" A ( M - 1 )(0 O)Tf)) O)Ti |0 1)T()) 1)T2 |0 2)Ti) 2)T2 ). (77) 

We first calculate the maximum eigenvalue A^T 2 ) and the corresponding left and 
right eigenvectors, \Z L ( M = 2 ) (p- TQ ; a Tl , cr T2 , a TZ ) and \/ R ( M=2 ) (er To , a Tl , cr T2 , 0> 3 ) , by solving 



(74) 
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eigenvalue equations (J5S|) and (J57J) . We then calculate the density matrix p for the 
restriction of the basis states of T A ^ M=1+1 ^ and T B{M=l+l ^: 

p(a 

= E ^ M = 2 )K,a Tl ,a r2 ,a T3 )^ M = 2 )K,<,<,a r3 ) (78) 



oy ,cr T3 



= E V L(M=2) (a T0 , a Tl , a r2 , a T3 ) ^ R ( M = 2 ) (a ro , a Tl , < , < 3 ) . (79) 

The left and right eigenvectors v L (o~ Tl ,a T2 ) and v^(a Tl ,a T2 ) of the density matrix 
p{cr T1 , cr T2 \a' , a' ) are determined by the eigenvalue equations (J59j) and ([BO]). 
j-A(Af=i+i/2) j n new basis states is obtained by using the eigenvectors of the m 
largest eigenvalues of the density matrix, v B and v^ 2 as 

T A ^ M ~ l+1 l 2 >{(Jo jTO , O!0 |T1 ,2, °~Q,t-a\°~1,t i Cl,T3 l a 2,ri i2 ) 

= E E E ^" A( ' M ~ 1 ' > ( C '"0,to, O'O.ti |g"l,T , C"l,r 2 |c r 2,ri) g"2,r 2 ) 

fO.Tj ,O"0,t 2 ""2,^ ,CT2,t 2 °'l,7-2 

x 6(a , r2 a , r3 |ai ir2 ai i7 - 3 )t;f o ((To, ri ,ao, r2 )^ 2 (cT2, ri ,cT2, r2 ) (80) 

where a and a 2 corresponds to the new basis states defined by the m largest eigenvalues 
of the density matrix. 7" B ( A/=1 + 1 / 2 ) j s similarly given by 

j-B(Af-l+l/2) 3 \ai jTl , 0- liT4 \a 2 , Tl , «2,t 2 ,3, CT2,r 4 ) 

= E E E ^ B(M=1 H CT o,75 ) ^o ) 7 3 kl,'ni ) cri,T4l cr 2,73,cr2,T4) 

ct 0,t 2 ,O"0,t 3 0"2,t 2 i°"2,r 3 0"1,t 2 

x 6(ai iri CTi ir2 |a2, ri a2, r2 )wf (ao, r2 ,ao, r3 )^ 2 (<T2, 7 - 2 ,cr2, r3 ). (81) 
The QTM of M = 3 is obtained by using T A ( M=1+1 / 2 ) and T B ^ M=1+1 ^ as 

y(M-3) ; a 0)T12 , (J 0jT3 , a ,T 4 ,5 \ a 2,T : «2,ri, 2 , CT2,t 3 , «2,r 4 ,5 ) 

= E T B{ - M 1+1 ^ 2 \a 0tTi5 |cr ljT3 , <7 liT0 |cr 2iT3 , a 2)T4 5 , cr 2j . 



X q-A(M 1+1/2)^^ 0^0,^,2, C"0,t- 3 1^1,^-0, Cn,T- 3 |«2,^ 1 , 2 )- ( 82 ) 

The density matrices needed for the calculation of T A ( M = 2 ) and T B ( M = 2 ) are calculated 
from the eigenvectors of the QTM of M = 3: 

p(a 

71,2 > °~T3 l a ri i2 5 ^73) 

= 5- V L ( M =*\a T0 ,a T1 , 2 ,a T3 ,a Ti ,)V R ( M = 3 \^ (83) 



CTt "0' Qt 4,5 



p((7 

= E V^ M=3 )(a T0 ,a T1 , 2 ^ (84) 

jvi(Ar=2) anCl q-B{M=2) are obtained by using the eigenvectors of the m largest eigenvalues 
of the density matrix: 

^ M (o~0,T : ^0,71,2,3 I ^l.TD) <T l,T4l a 2,ri, 2 , 3 ) cr 2,r 4 ) 
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= Y Y Y^ = 2) ( a O,r ,OCQ jTh2 ,ao,r 3 \o-l,r ,0-l, T3 \a 2 , Th2 ) 

q 0,t 1 2'°"0, T 3 a 2.T 1 2'°" 2 > T 3 "'iiTS 

X 0(<Ti l7S <7i, T4 |<72,^<T 2lT4 )T;f ) (ao,Ti ia , ^rg)^ («2,ri, 2 , ^rg) (85) 



rj-B(M=2) l I I \ 

i l^CT 0jT2 , ao,T 3 ,4,5 I <T 1,T2 5 °1,7B I a 2,T 3 ,4,5 5 Cr 2,T 6 J 

= 51 51 13 TB(M=1+1/2) ( a 0,T4 l6 l O 'l,r 3; O-l,r 6 |o-2,r 3 ,a2,r4, 5) O-2,r 6 ) 

fO,T3,C«0,T4 5 T2,T3,a2,T4 5 CT 1,T3 

X 6 (cr , T 2 (T 0,r 3 kl,7a CTl.TS ) W f (°" 0^3 > a 0,r 4 , 5 ) Va 2 (O- 2,r 3 , «2,T4, 5 ) • (86) 

<j-b(m=2) sriown j n fig ure an d equation (|HH)l is obtained by shifting the imaginary 
time indices Repeating the above procedures, we extend the QTM in the 

/3-direction, and the QTM of large M is obtained within a restricted number of basis 
states. 

3.5. Thermodynamic quantities 

The iterative extension of the QTM in the /3-direction corresponds to a successive 
decrease in temperature, which is related to the Trotter number as T = 1//3 = l/(/3 M). 
At each temperature, the free energy of the system par site is obtained from the 
maximum eigenvalue A max . 

F(T) = -\n\ max /(2(3). (87) 

The entropy S is obtained by taking derivative of the free energy with respect to 
temperature 

S(T) = --^F. (88) 
The specific heat C is similarly obtained by taking second derivative of the free energy 

C(T) = -T-^F. (89) 

The magnetic susceptibility \s and the charge susceptibility Xc are also calculated from 
the free energy at slightly different magnetic fields h z and chemical potentials \x 

X, = - (90) 

X, = - %jF. (91) 

The local quantities such as the magnetization and the density of electrons are 
obtained by using the eigenvectors V R and V L of X m ax- The magnetization (Sq) and 
the density of electrons (n ) are given by 

(5S> = Tr{e-e E S5}/Z = (V L \S*\V R ) (92) 



\r 

(n ) = Tr{e^ H n }/Z = (V L \n , TO \V R ) (93) 
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Figure 12. Charge susceptibility of the one-dimensional Kondo lattice model with 
hole doping 5 in the insulating state at half-fiUing^T]. The solid circle on the vertical 
axis represent the results obtained by the original zero-temperature DMRG method. 
A s is the spin gap at 5 = 0. J is the antifcrromagnetic coupling between the conduction 
electrons and localized spins, t is the hopping integral. 



where Sq tq is the operator at site % = and imaginary time Tq. The magnetic 
susceptibility Xs is then obtained by applying a small magnetic field h z 

Xs = ton(S5)/h*. (94) 

Similarly, the charge susceptibility Xc is obtained by changing chemical potential \x 

Xc = lim (K) /t +A /t - (no) J/ Ap. (95) 

A/i— +U 

As a typical example, the charge susceptibility of the one- dimensional Kondo lattice 
model calculated by the present method is shown in figure 

The nearest-neighbor and next-nearest neighbor correlations are also obtained by 
using the eigenvectors V R and V L and the QTM T 

(S*S*) = Tt{B-P*S' Q Sl}/Z = iy L \Sl T JSljV R )/\ max (96) 
(S Z S Z 2 ) = Tr{e-? H S*S Z 2 }/Z = (V L \S* iT0 TS* 2jT0 \V R )/\ max . (97) 

In equation (f9T)J) we have used S* commute with the local Hamiltonian h^i+i. The 
internal energy of the system is calculated from the above expectation values and we 
can calculate the specific heat C by taking derivative of the internal energy E with 
respect to temperature 

C(T) = (98) 

In order to confirm the accuracy of the results, we need to increase the number of 
basis states m kept in each block and the Trotter number M of the QTM. The error 
caused by finite M is usually scaled by M -2 , and the correct value in the limit of 
M — > oo is obtained by the linear extrapolation with respect to M~ 2 . The error caused 
by finite m does not follow a simple scaling form in general. Nevertheless, the truncation 
error in the norm of the eigenvectors (V L \V R ) is given by 1 — Y^k=\ w k , where w k is the 
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Figure 13. Schematic diagram of the calculation of imaginary-time correlation 
functions with the eigenvectors (V L \ and \V R ). 

kth largest eigenvalue of the density matrix. The correction for the thermodynamic 
quantities will be estimated from the truncation error. 



3. 6. Dynamic correlation function 

The Green's functions and dynamic correlation functions are obtained from the 
imaginary time Green's functions and correlation functions through analytical 
continuation to the real frequency axis^2j- The imaginary time Green's function is 
calculated from the eigenvectors of the QTM as 

G(r 2j ) = -Tr{e-? H c irT (r 2j ) 4(0) }/Z 

= - (V L \ Cia (r 2j ) cl(0)\V R ) (99) 

where (V L \ and \ V R ) are eigenvectors of the largest eigenvalue X ma x, and t 2 j = j((3/M). 
The local dynamic correlation function XAs( r 2j) cm the /3-axis is also obtained by 

Xab{t 23 ) = Tr{e-? H Mr 2j ) B^j/Z 

= (V L \A(T 2j ) Bi(0)\V R ). (100) 

The frequency dependence of the imaginary time Green's function and correlation 
functions are obtained by Fourier transformation along the imaginary axis 

3 M 

%) 4E^G(t 2j ) (101) 

1V1 3=0 

3 M 

XAB^n) = E ^ nT2j XA B (r 2j ) (102) 

M 3=0 

where uj n is the Matsubara frequency that is n(2n + l)/3 for fermionic operators and 
lirn/ 3 for bosonic operators. 

The real frequency Green's functions and dynamic susceptibility are obtained by 
the analytical continuation to the real frequency axis. The Pade approximations or the 
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Figure 14. Quasi-particle density of states calculated in the one-dimensional Kondo 
lattice model with hole doping 5 in the insulating state at half-filling 17,. 



maximum entropy method EH are used for this calculation. The Pade approximation 
is a fitting of G{iu) n ) or Xab^u) by rational functions of frequency iui n , which are 
analytically continued to the real axis by iu Tl — > uj + i5. 

The maximum entropy method is based on the spectral representations 

/OO p — TUJ 

p( W )-— (103) 
-oo 1+6 K 

/oo i e _ ™ 

-ImxAB(u)- Z&; duJ ( 104 ) 
-oo 7r l — e 1 

where p(cu) = — ^^G(cu + i5) is the density of state. This method finds the best p(uj) 
and xas(^) that reproduce G{r) and Xab(j)- 

The dynamic structure factor Sab^) is obtained from the imaginary part of Xab(^) 
through the fluctuation dissipation theorem 

lm X AB(oo) = tt(1 - e-^)S AB {u). (105) 
3.7. Finite-T DMRG with fixed Trotter number 

Reliable analytical continuation is only possible when we have accurate imaginary time 
correlation functions. In order to obtain correct dynamic properties, it is inevitable to 
use the finite system algorithm of the DMRG and to refine the eigenvectors V R and V L 
of the QTM. In the finite system algorithm we reconstruct the blocks T A and T B by 
keeping the Trotter number fixed. We first extend T A until T B is reduced down to the 
initial block of 7" B ( M = 1 ) ) and then extend T B to refine its basis states. We repeat such 
sweeps until the eigenvectors V R and V L are converged. Figure El shows the quasi- 
particle density of states calculated from the eigenvectors obtained by the finite system 
algorithm of the DMRG. In this figure we find the closing of the gap at the Fermi level 
to = with the hole doping 5 in the Kondo insulator |17j. 
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3. 8. Conserved quantities of the quantum transfer matrix 

The eigenvalues of the non-Hermitian matrix are not restricted to a real number. This 
fact sometimes causes uncontrolled instability generated from small numerical errors. In 
order to suppress such numerical instability and to keep the accuracy of the calculation, 
a precise diagonalization of the non-Hermitian matrix is needed. For this purpose, it 
is efficient to reduce the Hilbert space by using the conserved quantities of the QTM, 
which are related to the conservation law in the system. For example, S* + S? +1 in 
the Heisenberg model is conserved at any imaginary time evolution generated by h ii+ i. 
This means the following equation is satisfied for each imaginary time evolution from 

T~2k+i to T2k+i+l- 

Si,r 2k+i + Si+l,r 2k+i = Si,r 2k+i+1 + Si+l,T 2k+i+1 - {106) 

This equation is translated to 

qZ _ qZ _ _qZ , qZ ( '\ 07] 

%T2k+i h^k+i+l i+l,T 2S ._H ' w i+l,r 2fc+i+ i V iw V 

and the following equation is obtained after the summation over k 

2M-1 2M-1 

E = E (los) 

3=0 3=0 

This result clearly shows that Qi defined by 

2M-1 

Qi = E ( 109 ) 

3=0 

is a conserved quantity of the QTM [51], lo^]. 

Similar calculations are possible for charge degrees of freedom. For example, 

2M-X 

Pi=Y, (HO) 

3=0 

is a conserved quantity of the QTM, when h iji+1 preserves number of electrons rii + n i+ i. 

By using these conserved quantities, we can classify the basis states of the QTM. 
Since there is no matrix element between different subspaces classified by Qi and Pj, 
the QTM becomes block diagonal. The maximum eigenvalue X max of the QTM is then 
determined by calculating the largest eigenvalue in each subblock. The X ma x of the QTM 
is obtained in the subblock of Qi = Pi = in most cases. The density matrix calculated 
from the eigenvectors of the QTM is also block diagonal, and we need only diagonalize 
the density matrix within each subblock. The reduction of matrix dimension improves 
the accuracy of the diagonalization and stabilizes the calculation. 



3. 9. Complex eigenvalues of density matrix 

We finally comment on the problem of complex eigenvalues. This is caused by the 
non- Hermit icity of the density matrix, which allows appearance of complex eigenvalues. 
These are usually generated by the artificial truncation of basis states, and the imaginary 
parts of the eigenvalues are typically much smaller than their real parts. The complex 
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eigenvalues always appear with their complex conjugate, and their eigenvectors are 
mutually complex conjugate of the other one. Thus we can define two real vectors, which 
span the Hilbert space generated by the eigenvectors of the two complex eigenvalues, 
from the real and the imaginary parts of the eigenvectors, provided that the imaginary 
parts of the eigenvalues are sufficiently small compared with its real part. After checking 
the dual orthogonality of the left and right basis states = Sy, we can continue the 
DMRG calculation using the basis states defined by the real vectors. 

4. Application to two-dimensional systems 

Electrons in two-dimensional systems exhibit various interesting phenomena. Although 
their properties have been extensively studied for a long time, many unresolved questions 
are still remaining. In order to study the ground-state and thermodynamic properties of 
the system, the quantum Monte Carlo method and the exact diagonalizations have been 
used. However, when the electron density is away from particle hole symmetric point, the 
negative sign problem arises in quantum Monte Carlo simulations. Although the exact 
diagonalization method provides rigorous results, the number of electrons in the system 
is limited to small numbers. Thus, we need a new reliable method for two-dimensional 
electron systems away from half-filling, and various applications of the DMRG to two- 
dimensional systems have been considered. Most of the applications use mappings on 
to effective one-dimensional systems. However, the mapping is not uniquely determined 
and many long-range interactions appear, depending on the system. We thus need to 
find appropriate mapping specific to each model. 

4-1. Mapping to the one- dimensional lattice model 

In this section we consider two-dimensional electrons in a magnetic field. In two- 
dimensional systems under a perpendicular magnetic field, the kinetic energy of the 
electrons is completely quenched and the remaining macroscopic degeneracy is lifted by 
Coulomb interaction. Such an interesting system is realized in quantum Hall systems, 
where various electronic states have been observed. Here we describe how the DMRG 
method is applied to quantum Hall systems [34*| IHTj. 

In order to represent the Hamiltonian in matrix form, we first define one-particle 
states ^nx{x,v) of two-dimensional electrons. We use the eigenstates of free electrons 
in perpendicular magnetic field and represent the wavefunction in the Landau gauge. 
The one-particle states are represented by the following wavefunctions and each state is 
uniquely specified by only two numbers N and X: 



where N is the Landau level index and X is the x-component of the guiding center 
coordinates of an electron in cyclotron motion. are Hermite polynomials and Cn is 
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the normalization constant. £ = {Tic/ eH) 1 / 2 is the magnetic length, which is the length- 
scale of the system. The guiding center X is related to the momentum k y through 
X = k y £ 2 . Since k y is discretized under the periodic boundary conditions, the guiding 
center X takes only discrete values 

X n = 27i£ 2 n/L y (112) 

where L y is the length of the unit cell in the y-direction. After we have fixed the Landau 
level index, each one-particle state is uniquely specified by the guiding center X n . Thus 
the Hamiltonian is mapped on to an effective one- dimensional lattice model. 

The ground state of two-dimensional electrons in a perpendicular magnetic field is 
determined only by the Coulomb interaction 

V(r) = —. (113) 
er 

The Coulomb interaction generates correlations between the electrons and stabilize 
various electronic states depending on the filling v of Landau levels. When the magnetic 
field is strong enough so that the Landau level splitting is sufficiently large compared 
with the typical Coulomb interaction e 2 /(e£), the electrons in fully occupied Landau 
levels are inert and the ground state is determined only by the electrons in the top most 
partially filled Landau level. 

The Hamiltonian is then written by 

H = S ^2 C li C n + -Z ^2 X! X! X! ^™i™2n 3 n4 C n 1 C n 2 c n 3 c n4 (H4) 
n " m ri2 «3 rt4 

where we have imposed periodic boundary conditions in both x- and ^/-directions, and 
5* is the classical Coulomb energy of Wigner crystal with a rectangular unit cell of 
L x x L^|33]. c\ is the creation operator of the electron represented by the wavefunction 
defined in equation (jlllj) with X = X n . A niri2n . irii are the matrix elements of the 
Coulomb interaction defined by 

A — S' 1 STS> ^ 

/i nirt2n37i4 ni+n2n3+n4 ^ ^ / , °m -n.A.q v J„jl1 



21)2 



L N (q 2 £ 2 /2) 2 



exp 



q^ 2 ■ / ,QxL 

— i(ni - n 3 



M 



(115) 



where Ln(x) are Laguerre polynomials with N being the Landau level index |HH]. 
8' n2 = 1 when ni = n 2 (mod M) with M being the number of one-particle states 
in the unit cell, which is given by the area of the unit cell 2irM£ 2 = L x L y . 

In order to obtain the ground-state wavefunction we apply the DMRG method 
described in section 2. We start from a small-size system consisting of only four one- 
particle states whose indices n are 1, 2, M— 1, and M, and we calculate the ground-state 
wavefunction. We then construct the left block containing local states of n — 1 and 2, 
and the right block containing n = M — 1 and M from the eigenvectors of the density 
matrices which are calculated from the ground-state wavefunction. We then add two 
local states n = 3 and M— 2 between the two blocks and repeat the above procedure until 
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Figure 15. Eigenvalues of the density matrix for two-dimensional system with 18 
electrons in 54 local orbitals. 



M local states are included in the system. We then apply the finite system algorithm of 
the DMRG to refine the ground-state wavefunction. After we obtain the convergence, 
we calculate correlation functions. 

The ground-state pair correlation function g(r) in guiding center coordinates is 
defined by 

eV e J 

where Rj is the guiding center coordinate of the ith electron, and it is calculated from 
the following equation 

q 2 ^ 2 ., ,q x L a 

«q • r z(m — n% 



2 v M 



9{r)= N(N -1) ^ £ 6XP 

X 5'ni-n my L y /2^\4iA. C n 3 Cm\^) (H7) 

where ^ is the ground state and A^ e is the total number of electrons. 

The accuracy of the results depends on the distribution of eigenvalues of the density 
matrix. A typical example of the eigenvalues of the density matrix for system of 
M = 54 with 18 electrons is shown in figure ITHj which shows an exponential decrease of 
eigenvalues w a . In this case the accuracy of 10 -4 is obtained by keeping only 200 states 
in each block. 



4-2. Conserved quantities in magnetic field 

The Hamiltonian written in equation ()114j) conserves center of mass of electrons G = 
J2i(X n )i, whose origin is the conservation of total y-momentum J2i{k y )i = J2i(X n )i/£ 2 . 
All the eigenstates of the Hamiltonian and the basis states of the blocks are classified by 
the quantum numbers, N e and G. These conservation laws restrict the basis states of the 
ground state, and the complete basis states of the extended system are not generally 
constructed from the eigenstates of the density matrix calculated from the ground- 
state wavefunction which has fixed numbers N e and G. Although this unexpected 
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Figure 16. The lowest excitation gap at various v in the lowest Landau level. 
Relatively large excitation gap is obtained at fractional fillings v — n/(2n + 1). The 
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Figure 17. Pair correlation function g(r) at v — 1/3 in the lowest Landau level. The 
length is in units of t. 



truncation is corrected by reconstructing the blocks using the finite system algorithm 
of the DMRG, we need to include additional basis states in the extended blocks to 
accelerate convergence in the infinite system algorithm of the DMRG. 

4-3. Diverse ground states in a magnetic field 

Here we present diverse ground states obtained by the DMRG method applied to the 
quantum Hall systems. In the limit of a strong magnetic field, the electrons occupy only 
the lowest Landau level N = 0. In this limit, fractional quantum Hall effect (FQHE) has 
been observed at various fillings [^ZJ - The FQHE state is characterized by incompressible 
liquid with a finite excitation gap |58j. The presence of the FQHE is confirmed by the 
DMRG calculations, where a relatively large excitation gap is obtained at various fillings 
between v — 1/2 and 3/10 [3JJ as shown in figure We clearly find large excitation 
gaps at fractional fillings v = 1/3,2/5,3/7,4/9 and 5/11, which correspond to primary 
series of the FQHE at v = n/(2n + 1). The pair correlation function at v = 1/3 is 
presented in figure IT71 which shows a circularly symmetric liquid state consistent with 
the Laughlin state 
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Figure 18. Pair correlation functions g(r) in guiding center coordinates, (a) The 
Wigner crystal realized in an excited state at v = 1/6 in the lowest Landau level. 
The number of electrons in the unit cell N e is 12. (b) Two-electron bubble state at 
v = 8/27 in the third lowest Landau level. N e — 16. (c) Stripe state at v = 3/7 in the 
third lowest Landau level. iV e = 18. 



In the limit of low filling v — > 0, mean separation between the electrons becomes 
much longer than the typical length-scale of the one-particle wavefunction. In this limit 
the quantum fluctuations are not important and electrons behave as classical point 
particles, whose ground state is the Wigner crystal. The formation of the Wigner 
crystal is also confirmed by the DMRG calculations at low fillings as shown in figure ITSl 
(a). The v dependence of the low-energy spectrum shows that the first-order transition 
to the Wigner crystal occurs at v ~ 1/7. 

With decreasing magnetic field, electrons occupy higher Landau levels. In high 
Landau levels, the one-particle wavefunction extends over space, generating effective 
long-range exchange interactions between the electrons. The long-range interaction 
prefers CDW ground states and Hartree-Fock calculations predict various CDW states 
called stripe and bubble e 59|. These CDW states are confirmed by the DMRG 
calculations as shown in figures HH1 (b) and (c). Although the CDW structures are 
similar to those obtained in the Hartree-Fock calculations, the ground-state energy and 
the phase diagram is significantly different |34j . The DMRG results are consistent with 
recent experiments[60j, and the discrepancy is due to the quantum fluctuations neglected 
in the Hartree-Fock calculations. 

5. Summary 

In this topical review, we have reviewed the DMRG method and its applications to 
finite temperature and two-dimensional systems. We first briefly reviewed the original 
DMRG method, which enables us to calculate the ground-state wavefunctions and low- 
energy excitations of large-size systems with controlled accuracy. The essential idea of 
the DMRG is the restriction of basis states using the eigenvectors of the density matrix 
calculated from the ground-state wavefunction. The truncation error of the wavefunction 
is estimated from the eigenvalues of the density matrix and we can extend the size of 
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system with controlled accuracy. 

This idea of the DMRG was applied to the QTM in one-dimensional systems, 
and a reliable method for the calculation of thermodynamic quantities at finite 
temperatures has been formulated. In this way the thermodynamic quantities are 
directly obtained from the maximum eigenvalue of the QTM, and temperature 
dependence is systematically obtained by extending the QTM in the /3-direction. The 
dynamic correlation functions are calculated from the eigenvectors of the QTM with the 
analytical continuation to the real frequency axis. 

We have also shown that the DMRG method is reliably applied to two-dimensional 
quantum systems in a magnetic field. We use the eigenstates of free electrons in the 
Landau gauge as initial basis states of two-dimensional electrons, and map the system 
on to an effective one-dimensional lattice model with long-range interactions. It has 
been shown that this method is successfully applied to quantum Hall systems and the 
existence of various ground states has been confirmed. 
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